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DFTS ON IRREGULAR GRIDS: THE ANTERPOLATED DFT 

VAN EM DEN HENSON * 



Abstract. In many instances the discrete Fourier transform (DFT) is desired for a data set 
that occurs on an irregular grid. Commonly the data are interpolated to a regular grid, and a fast 
Fourier transform ( F FT) is then applied. A drawback to this approach is that typically the data have 
unknown smoothness properties, so that the error in the interpolation is unknown. 

An alternative method is presented, based upon multilevel integration techniques introduced by 
A. Brandt. In this approach, the kernel, is interpolated to the irregular grid , rather than 

interpolating the data to the regular grid. This may be accomplished by pre-inultiplying the data by 
the adjoint of the interpolation matrix (a process dubbed anterpolation), producing a new’ regular- 
grid function, and then applying a standard FFT to the new function. Since the kernel is C°° the 
operation may be carried out to any preselected accuracy. 

A simple optimization problem can be solved to select the problem parameters in an efficient way. 
If the requirements of accuracy are not strict, or if a small bandwidth is of interest, the method can 
be used in place of an FFT even when the data are regularly spaced. 

1. The formal DFT and the ADFT. The DFT is defined as an operation 

that maps a length- A complex- valued sequence {j*o. iy_j} to another length- 

A' complex- valued sequence {j 0 . Ji AY-i} by the rule 

.Y-l 

( 1 ) si, = ^ * J * , for k - 0. I V - I . 

j=0 

As defined in (1), the DFT is performed on data that are presumed to be given 
on a regular grid, with constant spacing between the data points. Furthermore, the 

transform values {tq.si are also presumed to lie on a regular grid in the 

frequency domain. In many applications, however, the data for a problem are not 
spaced regularly. It is of some interest, then, to determine how a discrete Fourier 
transform may be computed for such a data set. To perform this computation, we 
develop and implement in one dimension an algorithm based on multilevel integration 
techniques outlined by Aclii Brandt ([2], [l]). The method presented here can also be 
developed for higher-dimensional problems. One application of this technique [G] is in 
the reconstruction of images from projections (inverting the Radon transform). 

To begin, it is necessary to decide what is meant by a Discrete Fourier Transform 
for irregularly spaced data. Therefore, the concept of a formal DFT is introduced, 
which is defined as follows: 

Consider any set of A ordered points in the interval [0, A"), satisfying 

0 < Jo < JT • • • < < V 

and suppose a vector- valued function (grid function) u(xj) is specified. The formal 
DFT is defined as the M - \I { + M 2 + 1 quantities 

9 -[ 

(2) up,) = £ “(*>"“ = Z F- --Ui < l < • 

j = 0 A 

where / is an integer, and M\ and M 2 are positive integers specifying the range of 
frequencies of interest. The formal DFT may be thought of as an approximation to a 
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selection (— M\ < / < A/ 2 ) of the Fourier coefficients of u(x). In this view, u( x) should 
be regarded as an A'-periodic function known only at the grid points Xj. 

It is desired that this sum be calculated to a prescribed accuracy, say c||u||i, where 
||u||i is the discrete L\ norm ||u||i = N~ l |u(x ; )|. Note that any grid spacing 

is allowed for the Xj (in particular, the spacing needn’t be constant), that there is no 
relationship between the integers A/], A/ 2 , and iV, and that there is no requirement 
that these integers have any special value, (such as being powers of 2). Calculating the 
sum in (2) directly would have a computational cost of O(MN) operations. Instead 
of forming this sum directly, though, an approximation to it will be computed, using 
an FFT to accelerate the computation. 

The procedure begins with the definition of an auxiliary grid, fl^, covering the 
range of values [0,-Y). Let N, be an integer, whose value will be determined shortly, 
and let the grid spacing h be defined by h = X/N *. Then the auxiliary grid consists 
of the points y T — (r — 1 )/z , for r = 1,2,... A\. 

Suppose that the value of some function g(y T ) on the regular grid is to be 
interpolated to the gridpoints x 3 by Lagrangian interpolation. We will identify the 
interpolation by specifying the degree of the polynomial to be used. Thus, p-degree 
Lagrangian interpolation is computed using a polynomial of degree p or less. For each 
Xj , the p + 1 nearest neighbors on the grid must be located. Let these points be 

designated */*(;, 0 )* y*(j,i) V*(j,p)- These points should be chosen modulo A\ so that 

a point near the limits of the interval [0, A ) may have neighbors near the other end 
of the interval. (This is justified because, as will be seen shortly, the function to be 
interpolated for the formal D FT is A'-periodic.) For each x 3 , the p + 1 Lagrangian 
interpolation weights are computed by 



(3) 



p 

w n( x j) = n 



m = 0 

m ^ u 






and the interpolation of the function g to the gridpoints i 3 is given 



p 

9(*j) ~ W ’‘( r j)*/(V«U,n)) ■ 

71=0 

Letting g be the vector of function values g(yk) and g be the vector of interpolated 
values g(x 3 ), the interpolation may be written in matrix form 



9 = [Zyl 9 

where is the A r X A r . interpolation matrix mapping a function on to the 

gridpoints {x^}. The entries of this matrix are 



r T r 1 _ / w n (jj) if K(j,n) = m 

1 ~ \ 0 else . 

We are now ready to compute an approximation to equation (2). The strategy 
will be to interpolate, for each u/, values of the kernel e~ lu>lT J from the auxiliarv grid 
• That is, equation (2) is approximated by 

N-\ 

( 4 ) u(ui) - u ( x i) e (^h x j) 

;=0 
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where 



c{ut,Xj) = W n( X j) e . 

n=0 

Let the column vector of exponential values e~ l ^ iyk he designated e/, and the vector 
of interpolated kernel values e(ui,Xj) be denoted i\ Then this interpolation can be 
written 



h = p;if/ ■ 

Notice, however, that c/ is to be used as the I th row of the matrix giving the kernel 
of the summation in (4). To compute <\ as a row vector, the adjoint of the matrix 
equation equation is needed, namely 

(nf = (r y r,) T = (n) T [r y } T . 

Let us define the 47-vector u = u (. *;/), the N-vector u = u(xj ). and the matrices giving 
the kernels of equations (2) and (4) as IF and IF, respectively. Let the 4/ x A\ matrix 
whose I th row consists of c~ l '*' iyn for the AA points on be designated IF (it is useful 
to observe that this matrix consists of 4/ consecutive rows of the standard D FT kernel 
for a uniform grid with AA points). Then the formal DFT is approximated 



u — \ r u 

^ \\u 

= U-{I X y} T t • 

This notation can be simplified slightly by denoting the vector created by multiplying 
u by [J^] 7 . as u. Since the matrix [Tf] T is the adjoint of the Lagrangian interpolation 
matrix, the process of computing u = [If] T u has been dubbed (interpolation. Then 
the approximation to the formal DFT is 

(5) u ^ \\ ii , 



which we call the Anterpolated Discrete Fourier Transform (A DFT). 

The ADFT , as a matrix multiplication, requires 0(4/ A.) operations. In general, 
X* will exceed X , so as a matrix- vector multiplication, the ADFT has no advantage 
over (2). If, however, X* is selected appropriately, the approximation can be computed 
quite rapidly. Let 4/, = max {4/j, 4 / 2 }. Then if AA is selected such that X m > 24/,, 
and at the same time X m is a number for which an F FT module exists, then the fast 
Fourier transform can be applied to compute the DFT summation 



FFT{rl} = -Lfr i i r e-**T‘/X- for / = -^+l,-4: + 2,...2l . 

T=0 11 1 



Recalling that h — A'/AA, it may be seen that the DFT summation therefore yields 
(1 /AA)u((2t/)/A'). Multiplying by A A thus yields a set of values that includes, as a 
subset, all the desired values of 

Computing the ADFT , then, consists of two phases: 

1. u is computed from u by anterpolation : u = [Iy] T u. 

2. u is computed from u by a Fast Fourier Transform. 
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2. Operation count for the ADFT. The cost of computing the ADFT con- 
sists of the cost of computing the interpolation weights, the cost of computing the 
vector u — [1*] T u, and the cost of the F FT on AT points. 

Computing the (p + 1 )N interpolation weights, w n (x ; ), by the formula in (3) is 
the cost of the computing the numerator, since the regular spacing on ft^ means 
that the denominators of w n (xj) are independent of j. To compute the numerators, 

p 

the product JJ ( x 2 — is computed for each Xj, requiring 2 p + 1 operations. 

m=0 

Then the n th interpolation weight can be obtained by dividing by the product of 

p 

(x m — with the precomputed denominator n (^(nj) “ requiring 2 

m = 0 

m^n 

operations for each of the p + 1 weights associated with the point Xj. The calculation 
of the weights thus requires 0(N(p- j- 1)) operations. It is important, however, to note 
that the calculation of the weights is dependent only on the relationships between the 
gridpoints { y *.} and and is independent of the data set, u(x 3 ). This means that 

if a known set of gridpoints {x ; } and a standard auxiliary grid ft^y are to be used 
repeatedly, the interpolation weights w n (Xj) may be precomputed and stored, and 
needn't be included in the cost of the algorithm. This will be assumed to be the case. 

The matrix [Ty} T is AT X A" and the data vector u is \ X 1, so the computation of 
u = [ly\ r u would be O(A'AT) if performed as a matrix-vector multiplication. There 
is, however, a much more efficient method. The index table k(j, n) can be stored 
along with the interpolation weights. For each Xj and for each n . the value of k(j, n) 
is the index of the n th interpolation neighbor that is used to interpolate from ft^. 
to the gridpoint x 3 . The periodic nature of the kernel being interpolated means that 
the interpolation is always to a gridpoint Xj in the center of the set of p interpolation 
points (as is well known. [5], the Lagrangian interpolation is better behaved when this 
is the case). If p is odd, then x 3 always lies between y^- £-i_j and y ^ y while if p 

is even, then y K ( J% t) is the closest gridpoint on ft ^ to xj. 

Computing the vector u is then very easy, and may be done in 2A’(p+l ) operations, 
according to the algorithm: 



1. Initialize u(y 7 ) = 0 for all y T g ft^ . (1 < r < AT) 

2. For j = 0,1, A T - 1 

For n - 0, 1, . . . ,p 

Set “(yicfj.n)) — fi (»K0,n)) + W n (Xj)u(Xj). 



Having computed the values of u, consider now the cost of the F FT portion of 
the ADFT. This is simply the cost of an AT-point F FT . In the next section criteria 
for choosing AT will be determined. For now the only requirement is that an F FT can 
be computed on a vector of length AT. As such, AT must have factors for which F FT 
modules are available. For the purpose of an operation count, however, it is easiest 
to assume that AT is a pow’er of 2. Indeed, we shall see that w*e have great flexibility 
in our selection of AT, and since powers of 2 or 4 produce the most efficient FFT s, 
this is a good assumption. In this case, the cost of the FFT portion of the ADFT is 
0 ( AT log 2 AT). 
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The costs of the ADFT can now be computed. In terms of data storage it requires 
four arrays. One is the Appoint vector containing the input data, u . In addition, an 
Appoint complex vector is required for the input and output of the F FT. Assuming 
that the weights are precomputed and stored, two auxiliary arrays are necessary, an 
N X (p + 1) real (or double precision) array holding the interpolation weights w n (xj), 
and the N x (p + 1) integer array of indices, k(j, n). 

If the operations of multiplication and addition are counted equally, and if the 
weights and indices are pre-stored, the operation count is C\N „ log 2 A\ complex oper- 
ations for the F FT portion of the algorithm, where C\ depends on the choice of F FT 
algorithm. The computation of u entails 2A r (p+l) operations that are real or complex 
according to whether u is real or complex. Counting both phases of the algorithm, 
the operation count of the ADFT is 



C i A , log 2 A . + 2 A (p + 1 ) . 



This should be compared with the operation count of the formal D FT . which is 
0(M X). The computation of the ADFT is more efficient provided M is larger than 
2{p -f 1 ) + (AC log 2 A\)/A\ a condition that will generally occur in practice. 

3, Error analysis for the Anterpolated DFT. One of the attractive fea- 
tures of the ADFT is that the interpolation is performed on the kernel, which has 
known smoothness properties, rather than the data set. which generally has unknown 
smoothness properties. Since interpolation error depends on the smoothness of the 
interpolated function, the error committed by using the ADFT is relatively easy to 
analyse. 

Consider the error in /> degree Lagrangian polynomial interpolation, when the 
interpolation is from a set of + 1 gridpoints that are equally spaced. Let these 

gridpoints be designated sp- A function /(x), whose values are known at 

these gridpoints, is to be interpolated to the point x E [scosp]* Let x = £ 0 + th, where 
h is the gridspacing, and t E [0 ,p). The approximation to /(£ 0 + th) is the value of 
the Lagrangian interpolation polynomial P p (£ 0 + th), 

p 

Pp{x) = E w ‘(- r )/(s.) where w,(i) 

i = 0 

Defining 7r p (t) = t(t - l)(t — 2) . . . (t - /;), and ( E [<fo - sp] - the error in the interpolation 
is bounded by 



TT ^ 1 ~ ^ 

II (c _ c \ * 

m = 0 
m i 



( 6 ) 



l/(G> + th) — P p( so + th ) | < 



\- p (t)\\h^\Q p+i 

(P+1)! 



where Q P +\ = max^ 0 ^ ] |/^ P+1 HC)I* See [8], pages 264-270, for a derivation of this 
error term. 

It is useful to bound this error more precisely. To do this, we examine the behavior 
of the factorial polynomial 7T p (t). This polynomial has been well-studied, and many 
results can be found in various numerical analysis texts, ([5], [8], [9], [11]). These 
results, however, are developed for the case that x can be anywhere in [£co£p]- In the 
present case the interpolation is always to the center subinterval. Thus for p odd, 
t 6 while for p even, either t E [§, § + 1] or t E [§ — 1 , § ]. 
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To shorten the discussion, assume that p is odd. This is the most common case, 
where p = 1 gives linear interpolation and p = 3 give cubic interpolation. Similar 
results can be obtained for p even. Consider the following lemma, the proof of which 
may be found in [6]. 

Lemma 1. If p is a positive odd integer , then 

MOI „ i 

max — - — — < — — r . 

£6[ *y!.£$i](j>+l)! 2 p+1 

This result can be used to find an error bound for the ADFT. In this case the 
functions being interpolated are 

V V V 

e-’u"* for / = + 1 ,-— + 2 ,...— . 

2 2 2 

From this set of functions, the only ones whose values are of interest are those for 
/ between -M\ and A/ 2 . Recalling the definition A/* = max { M \ , A/2}, the largest 
absolute value among the frequencies of interest is t j\/ = ( 2 rr A /, ) / X . Therefore 



max 



0 *+ l ( 

v 






x ) 



P.W 



.P+1 



Inserting this and the bound from Lemma 1 into equation (G) gives 

(7) k —,- £( _ WI < 



which is then used to obtain 

I « I 1 1 ' u | = | £- V = “’ u(Xj) (- ‘~‘ T > - T-'jo u(Xj)t(~l.Tj) | 

< T.j= ~o 

< (^a// 2) p+1 E ; =o , I“(^)I 



= (/io £ ’ A/ /2)p +1 A'||u||, . 

Finally, substituting h = X/X m and = ( 2 xM m )/X establishes the desired error 
bound. The error in the ADFT approximation to the formal DFT is bounded, for 
-A/i < / < A/ 2 , by 



/ \f TT \ P+ 1 

(8) |u(^)/) - [U'fi](^;)| < A r ||u||i, 

where **;/ = 2 x 1 /X, and the ADFT ( 5 ) is computed using an FFT of length A r .. 
Since the bound holds for all desired values of ^7, it then follows that 



(9) 



« - [w'5)IL £ 
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where || • is dpfined as the maximum absolute value in the vector. It is also worth 
noting that an error bound for any desired frequency can be obtained by replacing 
c with u>i in the derivation, leading to 

(10) |a(wi) - [IV’u](uJ,)| < .V||u||i ■ 

This is especially useful information for those occasions when only the low frequency 
components are of interest, or when the accuracy required of the approximation is 
greater for the low frequencies than for the high frequencies. 

4. Selection of p and X m . The error bound just derived is useful in that it 
provides a way to select the operational parameters X m and p. Recall that the goal 
is to calculate an approximation to u to some prescribed accuracy, | u - U'u| < c| | u| | . 
In practice we will want to make the error small, so it will be assumed that c <C 1. 
Comparison with (8) gives tlm requirement 

(^r * f 

which may be written as 

(11) .V. > M.~ • 

For a given formal D FT. the values of M m . X . and ( are considered to be part of 
the problem specification. To ensure that the specified accuracy is obtained, it is only 
necessary to select integers X m and p so that (11) is satisfied. Naturally, there may be 
many combinations of parameter values that achieve this goal. The parameters should 
therefore be selected to fullfill some other desirable property as well. Specifically, they 
should be selected also to minimize the computational effort of the algorithm. 

To see how this may be accomplished, recall that the work involved in comput- 
ing the AD FT with X m points on the auxiliary grid and ;>degree Lagrangian 
interpolation is 



0(X m log 2 X m ) + 0(X(p + 1)) . 

The value of the constant on the 0(A r (p + 1)) term depends on whether the weights 
and indices are pre-stored, or calculated "on the fly*'. For the analysis that follows, 
we assume the weights and indices are pre-stored. in which case the constant is 2. 
The constant on the first term depends on several factors. F FT s generally have 
a complexity of {X/q)\og(X/q) for some number q > 1. If the data have certain 
symmetries, then a specialized F FT may be used for faster computation ([3]. [7], [10]). 
The variety of available F FT algorithms pursuades us to leave the constant on the 
first term as an unspecified parameter, C \ . 

The total work in computing the ADFT can therefore be written as a function of 
the two parameters X , and p. For a fixed problem size (A" and A/.), and a prescribed 
error tolerance c, the work in computing the ADFT to the required accuracy is 



(12) 



W(N.,P) = C i A’. log 2 N. + 2N(p+ 1 ) , 



and we seek an optimal parameters minimizing VP(iV,,p) over all combinations ( N mj p ) 
satisfying (11), if such a choice exists. 

Limiting cases may be determined by examining nearest neighbor interpolation 
(p — 0), as well as extremely high degrees of interpolation (p — * • oc). Substituting the 
limiting values of p into (11), and noting that equality will suffice to ensure that the 
required accuracy is attained, we obtain bounds for the selection of A r *, namely 



A/*7t < N m < 



A/.jt X 



e 



for all values of p > 0. 

The existence and uniqueness of optimal solutions are fairly easy to establish. 
\V(X m ,p) is continuous with respect to each of its variables, and both of the first 
partial derivatives are everywhere positive. This observation leads to 

Lemma 2. Let S be the set {{X.,p): X. > A/.jr(A7e) 1/(H ' I) }. and let OS be 

that portion of the boundary of S given by : X. = j. Then 

if {xo,yo) £ S, there exists a point (£,7/) E dS such that IF(£, 7/) < \V(x 0 ,yo). 

Proof: Since (x 0 .y 0 ) E S, the point (£, y 0 ) £ OS , where £ = M m x (7- ) yo + 1 . 
Furthermore, £ < j 0 . Then since the partial derivative of the work function with 
respect to X m is everywhere positive, \V(£.yo) < fF(^o<yo)- 

The utility of Lemma 2 is that the optimization problem can be rewritten as a 
problem in a single variable. Since for every point in S there is some point along OS 
that requires less work, it is only necessary to seek a minimum from the points of OS . 
This can be done by parameterizing X m and p as functions of a single variable. 



(13) 



*-(?)* 



Then on OS we find that 



(14) X m = M m 7:b and p+ 1 = log 6 

Since 0 < p < oc, the value of b is restricted to the interval (l,A T /c]. Substituting 
these expressions into (12), the work equation may be rewritten as a function of b 
alone 




(15) = C,.U.t6 ]og 2 (.\/.-6)+ 2.Vlog 6 (^-) , 

and the problem is to minimize (15) subject to the constraint 1 < b < (X/e). Once b 
is determined, the necessary values of X m and p can be obtained from (14). We may 
now establish 

THEOREM 1. There exists a unique value b that minimizes (15) subject to 1 < 
b leq(N/e). Therefore the work function 

\V(N m ,p) = Ci A r . log 2 X m + 2 X(p+ 1) 

has a unique minimum, subject to the constraints 

X m > M m 7r ^ ^ and 0 < p < 00 . 
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Proof: \V(b) is continuous and differentiable with respect to b on (l,iV/c], Dif- 
ferentiating equation ( 15 ) yields 

(Hi) ir'(6) = A’,ln(A' 2 6)-^|p . 

where 



A'. 



, A' 2 = 2 A/.jt , and A 3 = 2 .Y In (^-V 



For H*'(6) = 0 , then. 6 must satisfy 6(ln b ) 2 ln( A'2 6) = A'3/AV Now U 7 is also contin- 
uous and differentiable on ( 1 , A/c], and differentiating yields 



ir"(6) 




/ In 6 + 2 \ 
U 2 (In6) 3 J 



Since b > 1 we see that \V"(b) > 0 for all b 6 (l,A r /<], so any critical point in the 
interval must correspond to a local minimum. 

It is apparent that W'(b) cc as b — 1. Examination of the endpoint b = X/e 

reveals that since AT > tt, AN > and f < 1. we have that ir'(.Y/f) > 0. Further, 
IF"(&) > 0 implies that W'(b) has exactly one sign change in the interval ( 1 . .V/ c] . 
The point b 0 at which this occurs is therefore a global minimum for IF(6), and the 
value IE(AT,/j). where 



X. 



M.~ho and p 



lo S6 0 




- 1 



is the unique global minimum for W on OS. 

The values of AT and p obtained in this manner are real numbers. There is only a 
limited number of integers for which efficient F FT sexist, and Lagrangian interpolation 
requires p to be an integer. Further, this entire discussion has been predicated on the 
assumption that p is an odd integer, although a similar analysis can be made for p 
even. Once the theoretical values of AT and p are determined, they must bp modified 
to allow computation. There is some flexibility in this, but certainly selecting AT to 
be the first integer larger than M.nb for which an F FT exists, and choosing p to be 
the smallest odd integer greater than 



i°gb 




- 1 



will suffice. 

In order to find the optimal values of p and AT it is necessary to find the value of 
b satisfying 

(IT) 6(ln6) 2 ln(A’ 2 6)= . 

A 1 

While an analytic solution of this equation cannot be found, Newton’s iteration may 
be used. Table 1 displays optimal parameters AT and p for several combinations of 
AT A/,, and c. 
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32 


8 


.1 


48.7 


7.7 




128 


32 


.1 


193 


9.9 


32 


32 


.1 


142.6 


15.5 




128 


64 


.1 


325.7 


13.8 


32 


64 


.1 


257.6 


22.2 




128 


128 


.1 


570.7 


19.4 


32 


8 


.01 


52.9 


9.8 




128 


32 


.01 


206.7 


12.1 


32 


32 


.01 


150.1 


19.1 




128 


64 


.01 


344.1 


16.6 


32 


64 


.01 


267.8 


27.1 




128 


128 


.01 


595.6 


23.1 



Table 1. Optimal parameters N m and p computed for various problems. 

5. An AD FT Example. To illustrate the ADFT , consider the problem of com- 
puting the formal DFT of the function u(x) = [( 7r - x)/i r] 2 , sampled on an irregular 
grid. The irregular grid consists of N — 128 points x ; randomly spaced in the interval 
(0,27 t). Since the extent of the interval is 2 ti\ the frequencies u ;/ are just the integers 
/, and the formal DFT is 

A r -i 

(18) «(/) = Y. u{xj)e- ilx >, — G4 < / < 64 . 

j=o 

The sampled data are shown in the top of Figure 1. The real part of (18) is plotted 
on the bottom of Figure 1. The ADFT was used to approximate the formal DFT , 
with values of X m = 128, 236, 512, and 1024. Figure 2 displays, for each choice 
of A’,, the absolute value of the error | a( /) — [lFu](/)|, plotted as a function of /. 
Linear interpolation (p = 1) was used in each case. Note that increasing the value 
of X m produces a noticeable decrease in the error, and that the error increases with 
increasing wavenumber, as might be inferred from (10). Figure 3 displays the effect of 
using different values of p for fixed X m . It may be seen that the error decreases rapidly 

as p is increased. Equation (9) predicts that the error should decrease at least as fast 
/ w \ p+i 

as ( decreases as p or X m are increased. Table 2 gives both the infinity and 

Li norms of the error \u(l) - [U’u](/)| for several values of p and each of X m = 256 
and X. = 512. For .V. = 256, the error bound decreases by 0.6169 each time p is 
increased by 2. The experimental error is diminished by a factor of approximately 
0.3 as p is increased from 1 to 3, and by a factor of approximately 0.4 with each 
succeeding increase, better than the theory predicts. Similarly, for A T * = 512, the 
theoretical bound decreases by 0.15421 as p is increased by 2, while the experimental 
decrease is approximately 0.11 for each increase, a slightly better result. Numerical 
experiments on numerous other irregularly sampled functions, with various degrees 
of smoothness, produced similar results. In these experiments the ADFT behaved 
in a similar fashion as it did for the function discussed above. There is dramatic 
improvement with increasing values of A«, and p. As might be expected, the error 
diminished faster with smooth functions than discontinuous functions. 
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256 
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1.20663 


0.434861 




512 
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0.290501 


0.109943 


256 
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0.357889 


0.116080 




512 


3 


0.029351 


0.008315 


256 


5 


0.144478 


0.039136 




512 


5 


0.003360 


0.000817 


256 


*7 

i 


0.061305 


0.014983 




512 


7 


0.000419 


9.66385 e -05 



Table 2. Errors of the ADFT for various values of A * and p. 
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6. Some Open Questions about the ADFT. Like the continuous Fourier 
transform, the DFT has several important properties, such as linearity, the convolu- 
tion and correlation properties, the shifting property, the modulation property, and 
Parseval’s relation. To what extent these properties hold for the ADFT is an open 
question. The linearity holds can be established immediately, by noting that both the 
formal DFT and the ADFT can be written as matrix operations, so they are linear 
operators. Certain symmetry properties are easy to establish. For example, applying 
the ADFT to a real-valued vector will yield a conjugate symmetric result, that is 
u(cj) = u(— u>), because the vector [l*] 7 u is real-valued, and because the ADFT is 
computed by applying the F FT operator to this vector. The DFT , and therefore the 
F FT , maps a real vector to a conjugate symmetric vector [4]. Applying the DFT to 
data vectors w r ith other symmetries (even, odd, quarter- wave, etc.) yields output vec- 
tors with other types of symmetries [10]. It is natural to ask w'hich of these symmetry 
properties are inherited by the formal DFT or the ADFT . It seems reasonable to 
postulate that if the irregular gridpoints are symmetrically disposed and the function 
u(xj) is symmetric then the symmetry property of the DFT might be inherited by 
the formal DFT and the ADFT. 

Ail important question is: How is the formal DFT related to the continuous 
Fourier transform? That is, to what extent, and with what error, does the formal DFT 
approximate the FT1 Answering this question may prove to be a lengthy process. 
Many related questions will also arise. For example, how does the sampling theorem 
apply to an irregular grid? What frequencies can be represented accurately, and what 
constitutes aliasing? Is there some analog to the Poisson summation theorem? Many 
problems feature irregularly spaced data, so it may be assumed that these questions 
are of some interest. 
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Fig. 1 . The function f(z) = [(x - r)/ tj 2 sampled on an irregular grid, and its formai DFT. Only 
the read part of the formal DFT is plotted. 
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Fig. 2. The absolute vaJue of the error of the AD FT. plotted .as a function of wavenumber. Each 
graph corresponds to a different choice of .V., while linear interpolation {p = 1) is used for all. 
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Fig. 3. The absolute value of the error of the AD FT, plotted as a function of wavenumber. S m = 512 
for each graph, while several values of p are used. 
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